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Abstract: We present a new normalization condition for the axial current, derived from 
the PCAC relation with non-vanishing quark mass. This condition is expected to reduce 
mass effects in the chiral extrapolation of the results for the normalization factor Z\. 
The application to the two-flavor theory with improved Wilson fermions shows that this 
expectation is indeed fulfilled. Using the Schrodinger functional setup we calculate Zx{g$) 
as well as the vector current normalization factor Z\{g^) for (3 = 6/g$ > 5.2. 
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1. Introduction 

In recent years Wilson's formulation of lattice QCD [1] has matured to a stage where 
simulations with light dynamical fermions are within reach [2-4]. For light quarks some of 
the most interesting physics at low energies is associated with the pseudo-scalar sector of 
the theory, whose dynamics is governed by chiral symmetry. It is therefore important to 
study how and to which extent this symmetry can be realized in our chosen regularization. 

For Wilson fermions the (local) isovector axial current is not associated with a sym- 
metry of the lattice action and consequently it does not satisfy continuum Ward-Takahashi 
identities. However, the latter can be restored up to cutoff-effects (i.e. powers of the lattice 
spacing) through a finite rescaling of the axial current [5]. Its normalization factor Za(^q) 
is exactly obtained in this way by enforcing one particular continuum Ward identity at 
finite lattice spacing. 

Although a perturbative estimate of is available, most of today's simulations are 
performed at bare gauge couplings g$ ~ 1, where bare perturbation theory cannot be 
expected to work. To keep systematic effects in physical observables under control, it is 
therefore mandatory to determine non-perturbatively. 

The approach used by the ALPHA Collaboration in the quenched approximation [6] 
can obviously also be applied to the dynamical case. Here we improve this method by de- 
riving a normalization condition at finite quark mass, which provides a technical advantage 
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when extrapolating the results to the chiral limit. Using the Schrodinger functional setup 
we compute Za(<7q) f° r the relevant range of bare gauge couplings. In addition, we check 
the matching of our non-perturbative estimate with 1-loop perturbation theory through 
simulations at large f3. Systematic effects and in particular 0(a 2 ) uncertainties arising 
from a variation of the normalization condition are also considered. These turn out to be 
rather large for (3 < 5.4. 

The paper is organized as follows. After recalling the axial Ward identities we review 
the normalization condition used in [6] before introducing the new "massive" condition. 
Numerical results for Za and Zy are presented in Section ||, where we also provide inter- 
polating formulae. In Section || our conclusions are summarized and we discuss possible 
applications of the results. 

2. Theory 

We consider an isospin doublet of quarks with mass m and proceed formally in the contin- 
uum theory. The notation employed follows that of [7] and a pedagogical introduction can 
be found in [8]. 

2.1 The axial Ward identity in the continuum 

By performing local infinitesimal transformations of the quark and anti-quark fields in the 
Euclidean functional integral one derives the Ward identities associated with the flavor 
chiral symmetry of the action. An axial transformation gives the partially conserved axial 
current (PCAC) relation 

(d^A^O) = 2m(P a (x)0), where (2.1) 
= V^h/^lrXx) . pa ( x ) = #c)75^t<V(z) (2.2) 

are the axial current and the pseudoscalar density, respectively. Here r a are the Pauli 
matrices acting on the flavor indices of the quark fields. The relation (^]l|) holds for any 
operator O built from the basic fields in a region not containing the point x. 

Let 7Z be a space-time region with smooth boundary dlZ where the symmetry trans- 
formations are applied. The axial current A b v (y) is inserted as an internal operator (y£TZ). 
If O cxt denotes a polynomial in the basic fields outside this region the integrated form of 
the axial Ward identity is 

/ da^x)(A-(x)Al(y)0^) 
Jan N ' 

-2m [ d 4 x /p a (x)A b u (y)O cxt \ = ie abd (v u d (y)O cxt \ . (2.3) 

The right-hand side of fl2,3| ) originates from the variation of the internal operator A h v (y) 
since under chiral transformation the axial current linearly combines with the vector current 

Vuiv) =^{v)iu\rU{y) ■ (2.4) 
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In the relation (^^) we set v = and choose the region 1Z to be the space-time volume 
between the hyper-planes at yo~t and yo+t. After introducing a spatial integration over y 
and contracting the isospin indices with the totally antisymmetric tensor e abc we arrive at 



d 3 y f da^x)e ahc (Al{x)A b {y)O c 
Jan x 



-2m / d^y / d 4 x e abc 
K 



P a (x)A b (y)O cxt )=2i J d 3 y(v c (y)O cxt ) . (2.5) 



Note that the integral over 1Z includes a contact term at the point x = y. However, 
power counting and the operator product expansion tell us that the correlation function 
multiplying the mass has no non-integrable short-distance singularity. 

The above equation can be simplified by combining the two contributions from the 
surface integral over dlZ. We first notice that with periodic boundary conditions in the 
spatial directions J d 3 xdkf(x) = for any f(x) and thus the integrated form of the PCAC 
relation (^]l]), now written as an operator identity, reads 

Jd 3 xd A%(x) = J^xd^ix) = 2mjd 3 xP a {x) . (2.6) 

Integrating this relation from yo — t to yo then results in 

/r ryo r 

d 3 xA^(y -t,x) = d 3 xA a (y ,x) -2m / dx d 3 xP a (x ,x) . (2.7) 

With the region 1Z defined as above, the integration over the lower surface in (|2.5| ) involves 
current insertions at yo~t and yo, which can be shifted to yo an d yo + t by using (twice) 
the partial conservation of the axial current in the form of (^?j). The final result is 

Jd 3 yd 3 xe abc (A a (yo+t,x)A b (y)O cxt ) 
-2m y*d 3 yd 3 x J VO+ dxoe abc (p a (x)A b ) (y)O cxt )=i fd 3 y (v c (y)O cxt ) . (2.8) 



Note that this formal manipulation is not invalidated by the contact term in ( |2.5| ), which 
now appears at the lower integration limit. 

2.2 Normalization of the axial current on the lattice 

A normalization condition for the axial current on the lattice is derived by demanding that 
eq. ([^) in terms of the renormalized currents holds at non-zero lattice spacing. While 
this assumes the knowledge of the renormalization factor of the local vector current, it 
will become clear that with our choice of O ex t the latter can be easily computed. In the 
improved theory, eq. (2.8) then defines the normalization factor up to 0(a 2 ) uncertainties. 
We construct the relevant matrix elements in the framework of the Schrodinger functional 
[9, 10]. The notation is taken over from [7], to which we refer for any unexplained symbol. 
Our starting point is the improved axial current 

(A.r^A^ + acAd^, (2.9) 
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where denotes the symmetric lattice derivative. In a mass-independent renormalization 
scheme the renormalized (improved) current takes the form [7,11] 

(A R )« = Z A {1 + & A am q )(Ai)« . (2.10) 

Here m q = mo — m c is the bare subtracted quark mass, which on the lattice we need to 
distinguish from the current quark mass m, derived from a discretized version of ( |2.1[ ). 
In implementing fl2.8| ) on the lattice, we choose the isovector scalar 

O c ext = --^e cde O' d O e (2.11) 
as external operator. It is built from the zero momentum sources 

O a = « 6 Eu,vC(u)75^C(v) and 

/a = « 6 Eu,vC'(u)75^C'(v), (2.12) 

where £ and (' are the quark fields at the SF boundary xq = and xo = T, respectively. 
The isospin index of the external operator is contracted with the free index in ( |2.8| ) . As was 
shown in [6] , isospin symmetry in the form of a vector Ward identity implies that with this 
setup the right-hand side of eq. ( fO| ) simplifies to the boundary-to-boundary correlation 
function 

fi = ~^(0' a O a ) (2.13) 

up to 0(a 2 ). This gives the vector current normalization condition, which we will discuss 
in more detail in the next section. 

For the normalization condition derived in [6] ("old" condition) the mass m q was set 
to zero in ( [2.8; ) and ( |2.10| ). Since the normalization of the boundary fields cancels as they 
appear on both sides of ( [2.8D , with the external operator ( [2.11 ) the old normalization 
condition can be written as 

Z 2 A f{ A (y +t, yo) = fi + 0(a 2 ), (2.14) 
with the correlation function 

Aa^Vo) = --^^e afec e cde (0" i (A I )g(x)(A I ^(y)0 e ) . (2.15) 



In [6] the choice yo = t = T/3 was made in order to maximize the distance between the 
current insertions and thus reduce lattice artifacts. Since the mass term in the Ward 
identity is neglected entirely, an evaluation of this normalization condition at non-vanishing 
quark mass leads to errors of O(rom) in addition to the usual 0(a 2 ) errors. Here we have 
used ro, introduced in [12], as a typical low-energy scale. 

2.3 New normalization condition 

We now repeat the above steps using the full axial Ward identity to derive a normalization 
condition for the axial current. The additional term in eq. (|2.8|) proportional to the quark 
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mass results in a new correlation function /p A (2/o +t, 2/o)> which is similar to ( 2.15| ) but 
involves a temporal sum 

7 yo+t 

jp A (2/o+t,2/o) = E ^(^o)E ea6CeCrfe (° ,<ipa ( ;K )(^)o(y) oe ) » ( 2 - 16 ) 

zo=yo x,y 

where „(x ) = ( 1A X ° = y ° ° T V ° +t (2.17) 
I 1, otherwise , 

is introduced in order to implement the trapezoidal rule for discretizing integrals. Due 
to the aforementioned contact-terms on-shell 0(a) improvement does not apply to the 
expression ( |2.16| ). Therefore, although this new condition removes the O(rom) uncertainties 
in the determination of Z A at finite mass, other 0{am) ambiguities appear. In a quenched 
study [13] this residual mass dependence was shown to be very small. 

It follows from the PCAC relation that the product mP a renormalizes with the same 
factor as the axial current and hence the new massive normalization condition can be writ- 
ten as 

Z|(l+6 A am q ) 2 (/| A (yo+i,yo) - 2mfl A (y +t, y )) = fi + 0(om) + 0(a 2 ) . (2.18) 

In our calculation we neglect the coefficient 6 A , which is only known perturbatively. This 
choice only changes the 0(am) effects. The correlation functions for this new normalization 



condition can be evaluated with the same geometry as for (2.15), i.e. the volume integral 
covers the middle third of the temporal extension of the lattice. 

We now proceed to some technical aspects concerning the evaluation of f AA and / pA . 
Inserting the expression (|2,9j ) for the improved axial current into the correlation functions 
(|2l5|) and (|1|) we have 

Aa( x o, 2/o) = Jaa(xo, yo ) + oca 3o /pa(zo, 2/o ) + £>l /ap(^o, 2/0 ) 

+ a 2 ci^^/pp( 2 ;o,2/o), (2.19) 

yo+t 

fpA(yo+t,yo) =a^2 w(x ) /pa(3:o,2/o) +ac A ^/pp(x ,yo) , (2.20) 

£0=2/0 

where: / XY K Vo) = E e ahc e c ^(p" i X a {x)Y\y)0^ X,Yg{A ,P}. (2.21) 

Performing the Wick contractions for this correlation function, one finds that as a con- 
sequence of the isospin structure of the boundary composite fields O a and O la only six 
quark diagrams contribute. Those are shown in Figure |l[ Among them there are also two 
disconnected diagrams b) and c), where no propagator connects the xq = and x$ = T 
boundary fields. Exploiting the conservation of the axial current and making use of the 
operator product expansion, one can show [14] that in the continuum the diagrams b) 
and c) cancel for vanishing quark mass (the detailed argument is given in Appendix |A|). 
This implies that on the lattice they give only an O{o?) contribution to / AA and /p A in 
the improved theory. An alternative definition of Z A is thus obtained by using only the 
connected part of the correlation functions in fl2,18j ). 
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Figure 1: A graphical representation of the possible Wick contractions for the correlation functions 
fxY(xo,yo)- The gray diamonds indicate the insertions of Y and X at times yo and xo- 



2.4 The vector current 

Isospin symmetry, e.g. in the form of an integrated vector Ward identity, implies that on 
the lattice 

Z v (l + b Y am q )h(x ) = fi + 0(a 2 ) , (2.22) 
where fy{xo) is the Schrodinger functional correlation function 

M*o) = ^ E ^ abc (O' a V b (x)O c ) . (2.23) 

X 

Note that with spatially periodic boundary conditions, the term proportional to the im- 
provement coefficient cy [7,15] does not contribute to ( |2.23| ). In our implementation we 
also neglect by and extrapolate the simulation results to zero mass to obtain Zy from ( |2.22j ) 
in the chiral limit. The insertion point xq in ( 2.22j ) is taken to be T/2. The procedure 
is analogous to the one adopted in the quenched case [6]. There the slope (in am q ) of 
fi/(Zyfy) gives an estimate of by. This is no longer true in the unquenched theory as the 
improvement coefficient 6 g is different from zero [7] and thus the slope in am q is not due 
to 6y alone. 



3. Numerical Results 

We performed our simulations using non-perturbatively improved Wilson fermions [7, 16- 
18] and the plaquette gauge action. The clover coefficient c sw has been set to the value 
from [17, 18] and for the axial current improvement constant ca we have used the recently 
determined non-perturbative estimate [19]. Adopting the same setup as in the quenched 
computation [6,13] we choose T = 9/4L with periodic boundary conditions in space. The 
background field is set to zero. 

Concerning the algorithm, we employed the HMC with two pseudo-fermion fields as 
proposed in [20] and studied in the Schrodinger functional in [21] as well as the PHMC 
algorithm [22,23]. This choice has been motivated in [24], where PHMC has been shown to 
be more efficient at the coarsest lattice spacing, since there the spectrum of the improved 
Wilson-Dirac operator is affected by large cutoff effects in the form of unphysically small 
eigenvalues. 
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Summarizing the discussion from Section 2Ji , we choose the following non-perturbative 
definitions of the isovector normalization factors 



Z K {gl) = lim i ^[/I A (2T/3,r/3)-2m^ A (2r/3,T/3)l ^ , (3.1) 



m— >0 
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The quark mass m is defined as in [7] and we average it over three time-slices 1 around 
T/2. The same averaging is used to reduce the statistical error of Zy. Note that in the case 
of Za a similar time average would require additional components of the quark propagator. 
According to the discussion at the end of section we define another normalization factor 
Z A on through ( |3.1|) , where we now drop the disconnected quark diagrams in the correlation 
functions. 

As discussed in [6], we need to evaluate the normalization conditions on a line of 
constant physics, keeping all length scales fixed. This ensures that the 0(a 2 ) ambiguities 
in the normalization factors vanish smoothly when the perturbative regime is approached. 
In addition, the normalization conditions have to be set up at zero quark mass since we 
are aiming for a mass-independent renormalization scheme. 

To maintain a line of constant physics we need to know the three values of g$ that 
yield constant L for L/a = 8, 12, 16. Our value of L is determined by starting with the 
popular coupling = 6/g^ = 5.2 at L/a = 8. The readjustments of the bare coupling 
needed to achieve lattice spacings smaller by factors of 8/12 and 8/16 we take from the 
3-loop perturbative formula [25] 

= e -[ 9 o 2 -K)- 2 ]^o b 2 /( ^ )2] - bl ^ [l + ,[, 2 -(^) 2 ] + O((^) 4 )] , (3.3) 

q = 0.4529(1) , g < % . 

For the range of bare couplings involved this is also consistent with the non-perturbative 
coupling dependence of ro/a as checked in [19]. It should in any case be clear to the 
reader that small deviations from the constant physics line here influence the quality of 
improvement (size of remaining a 2 -effects) but do not imply any systematic errors of the 
continuum results. 

In addition to three matched lattice sizes L/a = 8, 12, 16 at /3 = 5.2, 5.5, 5.715, we 
simulated at three larger values of (3 and fixed L/a = 8, which corresponds to very small 
volumes. This was done in order to verify that our non-perturbative estimate smoothly 
connects to the perturbative predictions [26-28] 

Z A = 1 - 0.116458 gl + 0(g%) , (3.4) 
Zy = l- 0.129430 gl + O(g^) . (3.5) 

Our simulation parameters as well as the final results for Z\ and Zy are collected in 
Table [l], where we also include data from simulations at j3 = 5.29 and slightly mismatched 
volume. The latter are only used to qualitatively confirm the observed rapid change of 
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Table 1: Results for the chiral extrapolations of Z\ ( p.l[ ) and (3-2) and estimates for the 
critical hopping parameter k c . 



Za(^o) i n this region of the coupling. The results in Table [l| are obtained through an 
interpolation or slight extrapolation in the quark mass. The first error we quote for Z\ 
and Zy is statistical and the second represents our estimate of the systematic error, which 
originates from deviations from the constant physics condition. It is discussed in the next 
section. The complete set of simulations results and parameters is given in Table | m 
appendix [B|. 

As expected from the above arguments and already verified in the quenched case [13], 
for the new normalization condition the data exhibit very little dependence on the quark 
mass. Consequently, uncertainties in the location of the critical point do not propagate to 
the determination of Za- That this does not hold for the previously used normalization 
condition [6] can be inferred from Fig.g. There we show the chiral extrapolation at /3 = 5.2, 
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Figure 2: Comparison of the chiral extrapolation at j3 = 5.2 using the new and old normalization 
conditions for Za- 

where we had to use the PHMC algorithm in order to obtain reliable error estimates. With 
the HMC algorithm the error analysis would have been tainted by the rare occurrence of 

1 Four time-slices for odd T. 



-8- 




very small eigenvalues [24]. While for the new normalization condition the slope in am 
is consistent with zero, the estimate of from the old condition changes by 30% in the 
(small) mass range shown. We anyway see that for am < 0.02 all mass effects show a linear 
behavior. For j3 = 5.5 the extrapolation is similar to the one shown and at all other gauge 
couplings we can in fact interpolate from two simulations very close to the critical point. 

Our final results are shown in Fig. ||] as a function of g$. The errors plotted there 
are obtained by linearly summing the statistical and systematic errors. One can see that 
our data for both Za and Zy lie on smooth curves. These can be parameterized by the 
interpolating formulae 



= 1-0.918^+0.062^ +0.020^ 
1 0.8015 sg 



_ 1-0.6716,3 + 0-038M . (3 . 7) 
vyo; 1 - 0.5421 gl y ' 

To Zy we ascribe an absolute error of 0.005, whereas for Za the absolute error decreases 
from 0.01 at /3 = 5.2 to 0.005 at j3 = 5.7. The expressions in ( |3.6| , |3.7| ) have been obtained 
by a Pade fit, constrained by 1-loop perturbation theory ( |3.4| , |3.5| ). Though slightly mis- 
matched and hence excluded from the fit, also the (3 = 5.29 data (lighter color in Figs. 
|3| and |j) are reproduced by the interpolating formulae in both cases. For Z\ we find 
agreement at the 1% level with the results from [29], where isospin charge conservation 
is imposed in large volume for matrix elements of the local vector current among nucleon 
states in the bare coupling range (5 < 5.4. 

We note in passing that similarly to the quenched case [6] also here the use of mean-field 
improved perturbation theory [30] for ^a/v improves the 1-loop approximation. While the 
difference to our non-perturbative determination is small for g$ < 1, it rapidly increases 
in the range of physically relevant bare couplings. In particular, at a lattice spacing of 
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roughly 0.1 fm, corresponding to f3 = 5.2, our non-perturbative estimate of Za is almost 
20% smaller than the 1-loop value (10% for boosted perturbation theory). In the quenched 
case [6] this difference was roughly a factor two smaller at the same lattice spacing. 2 



3.1 Systematic effects 

Close to the continuum the dependence of the normalization factors on the lattice size is 
expected to be of order (a/L) 2 [6] in the improved theory. This implies that effects in Z\ 
and Zy due to deviations from the line of constant physics should be strongly suppressed. 

To check for these effects, at (3 = 5.5 and 5.715 the simulations closest to the critical 
point were repeated on smaller lattices (L/a = 8 at ft = 5.5 and L/a = 12 at ft = 5.715) 
in order to numerically assess the derivative of Zj^/y with respect to L. Moreover, we 
estimate the uncertainty in L (measured in units of L at /3 = 5.2) due to our approximate 
matching to increase linearly in ft up to at most 10% in the range 5.2 < ft < 5.715 (see [19] 
for a discussion of this estimate). We therefore assign a 6% error to L at ft = 5.5 and 10% 
at ft = 5.715. Together, this gives the systematic errors quoted in Table || through a linear 
propagation of the error. Our results confirm the expected small volume dependence. 

For the runs with ft > 7.2 the matched L/a would be extremely large. On the other 
hand the volume dependence should decrease like <7g(a/L) 2 as we are approaching the 
perturbative regime. In practice the simulations are thus performed at L/a = 8 and the 
systematic error is estimated from additional runs at the coarsest of these lattice spacings, 
i.e. at P = 7.2, by taking the difference of Z\ and Zy between L/a = 8 and 16. While for 
Za the volume dependence is hardly visible, in the case of Zy the larger volume {L/a = 16) 
results in a statistically significantly lower value. However, even this amounts only to a 
0.5% effect in the final result. In addition, we checked the dependence of Zy on the 
background field (BF), observing only a moderate effect when going from BF=0 to BF=A 
(used to define the running coupling g 2 in [31]) even for 5.2 < P < 5.4. 



3.2 Comparison with an alternative normalization condition 



As discussed in Section the disconnected diagrams in the correlation functions ( 2.21| ) 
are expected to contribute only 0(a 2 ) to Za- Hence they can be dropped to obtain an 
alternative normalization condition for the axial current. 

We carry out the previously discussed analysis also with this definition to obtain the 
results shown in Fig. ||, where they are compared to the interpolating formula (|3.6j ). The 
large difference for g 2 . > 1.1 has to be interpreted as a cutoff effect and we could indeed 
explicitly verify that it vanishes faster than linear in a. While this is consistent with 
the expectations from an improved theory, the magnitude of these lattice artifacts is still 
worrisome. These findings add further evidence that for a ~ 0.1 fm cutoff effects with 
dynamical (improved) Wilson fermions can be unexpectedly large [32]. 



2 In fact, in [6] Z£° u was considered, but the difference from our definition of Za was found to be negligible 
already for a ~ 0.1 fm. 
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Figure 4: Comparison of Z£° n and the interpolating formula (3.6) for Za- 



4. Conclusions and outlook 

In this work we have shown that in a lattice theory with two flavors of Wilson fermions 
normalization conditions can be imposed at the non-perturbative level such that isovector 
chiral symmetries are realized in the continuum limit. Since we are working with an 
improved theory, chiral Ward-Takahashi identities are then satisfied up to 0(a 2 ) at finite 
lattice spacing. 

The normalization condition was implemented in terms of correlation functions in the 
Schrodinger functional framework and evaluated on a line of constant physics in order to 
achieve a smooth disappearance of the 0(a 2 ) uncertainties. Through additional simulations 
at very small lattice spacings and volumes we verified that our non-perturbative definition 
approaches the perturbative prediction at small bare gauge coupling. 

Simulations were done at or near the critical point and owing to the new normalization 
condition, which keeps track of the mass term in the PCAC relation, any chiral extrapo- 
lations are extremely flat. Systematic effects due to deviations from the constant volume 
condition are also estimated and turn out to be small. The results are well described by an 
interpolating formula Za(<7 2 ) in the range of bare couplings considered. Enforcing isospin 
symmetry with the same programme, we obtain at the same time a non-perturbative de- 
termination of the normalization constant Zy{g§) of the vector current. Within about 1% 
it is in agreement with Zy determined in [29] but extends to weak couplings, where contact 
with perturbation theory is made. 

We found rather large 0(a 2 ) uncertainties in Za at < 5.4 by varying the definition 
of Za- Together with the algorithmic issues discussed in [24] these findings corroborate 
the worries expressed in [32] about the status of simulations with improved dynamical 
Wilson fermions at the currently accessible lattice spacings. This merely emphasizes that 
cutoff effects in physical quantities can be controlled only if a continuum extrapolation 
with several lattice resolutions is performed. 

The result obtained here is an essential step in the computation of the pseudo-scalar 
meson decay constant Fps needed to reliably convert the A parameter from [31] into phys- 
ical units. In the short term, together with data from [33] Za(<? 2 ) will be used in a fully 
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non-perturbative calculation of the strange quark mass following the strategy of [34] . 

Finally, the method employed here can also be used to obtain Z\ in the 0(a)-improved 
three flavor theory with either Iwasaki or plaquette gauge action [18,35,36]. 
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A. Contribution of disconnected quark diagrams 3 

In Fig. [l] the diagrams, which are related by an exchange yo *-> xq have the same isospin 
factors with opposite signs. This follows directly from ( |2.15 ), where such an exchange cor- 



responds to e a6c <-> e bac . Following [14], we now argue that in the massless continuum limit 
the contributions of the disconnected diagrams b) and c) to /aa(^o> Ho) cancel. Considering 
e.g. diagram b), its continuum version is proportional to 

/ A b {(x ,y ) = ja^ah{o' u Mo) d u{x)(Av) cs {y)O sc ) , (A.l) 

where we have introduced flavor indices u, d, s and c for the valence quarks, such that 
diagram b) is the only possible Wick contraction. Since the spatial insertion points x and 
y are integrated over, f^ A depends on xq and yo only. In fact, as we are in the chiral limit 
and the operators O and O' generate zero-momentum states, the axial current is conserved 
and hence the diagram is independent of the insertion points in the two regions xq < yo 
and xq > yo. If the two points meet, contact terms may arise, which we need to treat 
separately. 

To this end we restrict the spatial integration to |x — y| > e and let xo approach yo 
from either region. No contact terms can appear due to the finite spatial separation. In the 
limit e — > the contribution to f AA from the region |x — y| < e vanishes if the integrand has 
a divergence weaker than x — y[~ 3 . In this case we can safely take the limit and conclude 
that the order of xq and yo does not play any role. This would imply that the diagrams 
b) and c) have the same value and since their isospin factors have opposite signs, their 
contributions to /aa cancel. 



It is clear that the flavor assignment in (A.l) excludes a single quark bilinear as the 



leading contribution in the short distance expansion (for x — > y) of ^4o(^)^4o(y)- Hence, 



We thank Martin Liischer and Stefan Sint for their contributions to a clarification of this issue. 
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by power counting, the latter has (if any) a divergence weaker than |x — y|~ 3 in the limit 
|x — y| — > and the contribution from the excluded integration region vanishes. 

Since the correlation functions approach their continuum value with a rate proportional 
to a 2 , we can conclude that on the lattice the contribution from the disconnected diagrams 
is a cutoff effect of this order. 

B. List of simulation parameters and results 



In Table || the simulation results for Zx and Zy are collected. The number of measurements 
-^Vmeas explicitly contains the number of replica and r meas is the number of (unit length) 
trajectories between consecutive measurements. 



P ft L T N mcas ^meas 


am Z A Z%> n Zy 


5.200 0.13550 8 18 16-200 4 
5.200 0.13550 8 18 16-40 10 
5.200 0.13560 8 18 16-225 3 
5.200 0.13570 8 18 16-230 2 
5.200 0.13580 8 18 16-230 2 


0.01718(90) 0.7301(173) 0.8411(80) 0.7509(6) 
0.0159(11) 0.7186(295) 0.8455(108) 0.7497(14) 
0.01310(68) 0.7157(137) 0.8212(96) 0.7471(7) 
0.0088(11) 0.7134(126) 0.8302(70) 0.7447(8) 
0.00194(81) 0.7176(114) 0.8588(99) 0.7424(14) 


5.290 0.13625 8 18 16-50 2 
5.290 0.13641 8 18 16-120 2 


0.0031(18) 0.7527(102) 0.8103(167) 0.7507(19) 
-0.00512(61) 0.7540(124) 0.8378(73) 0.7490(12) 


5.500 0.13606 12 27 16-25 6 
5.500 0.13650 12 27 16-44 3 
5.500 0.13672 12 27 16-80 3 
5.500 0.13672 8 18 1-318 4 


0.02254(26) 0.8417(222) 0.8077(26) 0.7853(14) 
0.00758(27) 0.7987(153) 0.8100(45) 0.7738(8) 
0.00041(25) 0.7888(32) 0.8048(54) 0.7650(21) 
-0.00168(62) 0.8105(64) 0.8168(38) 0.7750(45) 


5.715 0.13665 16 36 1-106 2 
5.715 0.13670 16 36 1-54 2 
5.715 0.13670 12 27 4-62 2 


0.00194(57) 0.8142(135) 0.8079(31) 0.7827(11) 
-0.00060(69) 0.8004(26) 0.8120(30) 0.7793(20) 
-0.00100(34) 0.8021(38) 0.8182(18) 0.7861(18) 


7.200 0.13420 8 18 1-220 2 
7.200 0.13424 8 18 1-164 2 
7.200 0.13424 12 27 16-50 2 
7.200 0.13424 16 36 1-80 2 


0.00029(45) 0.8721(24) 0.8772(18) 0.8573(9) 
-0.00028(42) 0.8683(22) 0.8732(18) 0.8553(6) 
-0.00049(15) 0.8685(23) 0.8717(8) 0.8543(18) 
-0.00023(41) 0.8678(18) 0.8670(17) 0.8508(18) 


8.400 0.13258 8 18 4-40 2 
8.400 0.13262 8 18 4-45 2 


0.00023(40) 0.8990(28) 0.8956(16) 0.8839(15) 
-0.00183(42) 0.8998(25) 0.8953(13) 0.8826(7) 


9.600 0.13140 8 18 4-100 2 
9.600 0.13142 8 18 4-125 2 


0.00021(15) 0.9137(14) 0.9154(7) 0.9040(4) 
-0.00059(15) 0.9118(12) 0.9155(7) 0.9034(4) 



Table 2: Summary of simulation parameters and results for Za and Zy. 
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